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Abstract 



In this paper, we introduce a new nonlinear evolution partial diflterential equation for sparse deconvo- 
lution problems. The proposed PDE has the form of continuity equation that arises in various research 
areas, e.g. fluid dynamics and optimal transportation, and thus has some interesting physical and geomet- 
ric interpretations. The underlying optimization model that we consider is the standard £i minimization 
with linear equality constraints, i.e. min„{||M||i : Au = /} with A being an under-sampled convolution 
operator. We show that our PDE preserves the £i norm while lowering the residual \\Au — f\\2- More 
, importantly the solution of the PDE becomes sparser asymptotically, which is illustrated numerically. 

' Therefore, it can be treated as a natural and helpful plug-in to some algorithms for ii minimization prob- 

lems, e.g. Bregman iterative methods introduced for sparse reconstruction problems in Numerical 
experiments show great improvements in terms of both convergence speed and reconstruction quality. 

> 

O ■ 1 Introduction 

^ The sparse deconvolution problem is to estimate an unknown sparse signal which has been convolved with 

' some known kernel function and corrupted by noise. There are many applications of sparse deconvolution, 

' e.g. in seismology, nondestructive ultrasonic evaluation, image restoration, or intracardiac electrograms (see 

e.g. aaaal&a). 

There has been an intensive development of sparse deconvolution techniques and algorithms. Due to their 
. simplicity, least squares methods have been used to deconvolve sparse seismic signals 0, The major 

, drawbacks of least squares methods are their lack of robustness and the ill-poseness of the problem when the 

- number of unknowns is less than or equal to the number of equations. Many complementary methods have 

been presented, such as Tikhonov regularization 11 1 and the total least squares method [ij, [l^]- Another 



class of complementary methods to the classical least squares method is the £i-penalized models 

where £i norm of the unknown signal is used as an additional regularization term, allowing us to control the 

sparseness of the solution. 

In the recent burst of research in compressive sensing (CS) 0, EH following £i minimization 

problem (basis pursuit) is widely used (see e.g. [iqI [i|) 

mm {\\uh:Au^f}, (1) 
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where ||u||i :— \ui\ denotes the £i norm of u. In our paper, we shall adopt the above optimization model 
with A being a convolution matrix with some known continuous differentiable kernel function (e.g. the 
impulse response of the system or wavelet basis). When the kernel function is everywhere nonnegative and 
the sparse signal we want to recover is also known to be nonnegative, we can solve the constrained least 



squares problem min„{||Au — /||| : u > 0} as suggested by [20|. However, this model does not work for 
general kernel functions and signals with negative entries. Hence in this paper, we will stick to the above £i 
minimization model ([1]), although we believe our PDE can be also added to many nonnegative least squares 
approaches. 

To solve the constrained ogtimization ([T]) , standard second-order methods such as interior-point methods 
were introduced in [2l[ 22, [2J]. However, these methods, although accurate, become inefficient for large 



scale problems. To overcome this, large scale optimization methods have recently been developed to solve 
(P) or its siblings [l|, El S, 0, S, S El Q- These methods are very popular in CS, and all of 
them, especially Bregman iterations introduced to CS in have outstanding performance in solving CS 
problems, e.g. when the matrix A is taken to be a random submatrix of some orthonornial matrix, such as 
Fourier or DCT matrix. 

For deconvolution problems however, the above mentioned methods are usually less efficient than solving 
CS problems. The major reason is that each column vector of a convolution matrix A is highly coherent 
to the nearby columns, especially the ones that are right next to it. This makes the signal much harder to 
deconvolve if it has two spikes that are very close to each other. Therefore, in contrast to CS problems, 
the spatial information of u (i.e. the variable "x" of u{x)) is no longer irrelevant for sparse deconvolution 
problems, which motivates us to introduce what we call a spatial motion PDE into optimization algorithms 
to solve problem ([T]). 

Intuitively speaking, this PDE should spatially move spikes around such that the constraint Au = / is 
closer to being satisfied. However, the £1 norm of the solution should not be altered during the process. We 
will prove in Section [3] that our motion PDE indeed satisfies these properties. Furthermore, the PDE alone 
tends to sparsify the solution, which, together with its other properties, makes it an effective and powerful 
add-on to some of the optimization algorithms for CS problems mentioned above. In this paper, we shall 
combiiie^the motion PDE with Bregman iteration [l| (with the inner unconstrained problem ([3]) solved by 
FPC [30]) to enhance its performance. We note that it is also possible to combine the motion PDE with 



other iterative thresholding based algorithms (e.g. the linearized Bregman iterative method in 25| or the 



Bregmanized operator splitting method in 28[). 

The rest of the paper is organized as follows. We will first recall some basics of Bregman distance and 
Bregman iterations in Section [5] Then in Section [31 we introduce our motion PDE and explore its physical 
and geometrical meaning. In Section |4l we combine our PDE with Bregman iteration and propose our new 
algorithm. Numerical experiments are also conducted. Finally, we draw our conclusions in Section [SJ 

2 Background 

In this section, we briefly recall the concept of Bregman distance [s^ and Bregman iterations [l| . 
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2.1 Bregman Distance 

The Bregman distance [33|, based on the convex function J, between points u and f , is defined by 



D^j{u, v) — J(u) — J{v) — {p,u — v) 

where p € dJ{v) is an element in the subgradient of J at the point v. In general D^j{u, v) ^ D^{v, u) and the 
triangle inequality is not satisfied, so is not a distance in the usual sense. However it does measure 

the closeness between u and v in the sense that D^j{u, v) > and D^j{u, v) > D^{w, v) for all points w on the 
line segment connecting u and v. Moreover, if J is convex, Dj(u, v) > 0, if J is strictly convex D^j{u, w) > 
for u ^ V and if J is strongly convex, then there exists a constant a > such that 

D^{u, v) > a\\u — wjlj. 



2.2 Bregman Iterations 

To solve (IlJ Bregman iteration was proposed in Given — = 0, we define: 

yfe+i ^ argmin„e_R,. ~ A^lk'^lli - {u - u'',p^) + l\\Au - f\\l} 

pk+i ^pk _AT{Au''+^ - f). 

This can be written as, for J{u) = fi\\u\\i, 



(2) 



,fe+i 



argmin i?^ {u,u'') + -\\Au - fg 



As shown in see also 



34 



35|, the Bregman iteration © can be written as, for ={),vP = 0: 



arg min < + — /' 



k+l 

fk+l _ fk . r _ A^k+l 



k\\2 
2 



(3) 



r + f- Au' 



At each iteration, we will solve the problem ([3]) via the fixed point continuation (FP C) method propos ed in 
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36 



37 



38, 



32 



391 and 



[30j. Other ways of solving the unconstrained problem ([3]) can be found in e.g. 
the references therein. Our improvement will work with most of these solvers. Now altogether, we have the 
following FPC Bregman iterations solving problem ([T|): 



,k+l,N 



,k,l + l 



shrink(M''''+2 ^ ^) 



Here "shrink" is the soft thresholding function [4^ defined as 



(FPC) 



(4) 



shrink(x, ^) 



X — fi, if X > /i 
0, if — /i < 2: < 

X + i-i, if X < — /i. 
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Notice that, if we only perform 1 iteration for the FPC in i.e. N = 1, this gives us another algorithm, 
called Bregmanized operator splitting, that also solves problem ([1]) [^l- 



3 A Spatial Motion PDE 

In this section, we introduce and analyze a spatial motion PDE. The main purpose of this PDE is to spatially 
move and properly combine spikes in a sparse signal u, such that the residual \\Au — f\\2 is reduced while 
II M 111 is preserved. We will also present some interesting physical and geometric interpretations of the motion 
PDE. Throughout this section, u is understood as a function defined on the image domain $7, while A is a 
composition of an convolution operator and a spatial restriction operator with respect to a sampling region 
S, and is the conjugate operator of A. In other words, 

{Au){x)^Xs{x) I K{x-y)u{y)dy (5) 
Jn 

where K ]s a. differentiable convolution kernel and xs is the characteristic function of S. The other operators 
in this section should be understood as functional operators as well. 

3.1 Motivation 

To analyze the desired spatial motion of the signal, we start from a characteristic point of view, i.e. observing 
the moving trajectory of a point mass and describing its motion. Consider a simple sparse signal that consists 
of an isolated spike at location c: 

u = 5{x — c). 

We want to find the proper spatial motion that decreases the energy H{u) — ^||^u — /||2- To do this, we 
take the partial derivative of the energy with respect to the spatial variable c and obtain 

^H{u) = -V^{H'{u))l^^ = -VM^{Au - f ))l^^. 

The above identity means that if we want to minimize H{u), we can spatially move the spike along the 
direction \7 ^iA^ {Au ~ f))\,j.-^- In other words, the trajectory of the moving spike that minimizes H{u) can 
be described as u{x,t) = S{x — c{t)) where c{t) satisfies 

c'{t)^-VM'^{Au-f))l^^^^y 

Similarly, if the initial signal is a negative spike, i.e. u — —8(x — c), then its correct direction of spatial 
movement should be ^H{u) = V^iA^ {Au — /))|^^^. 

Based on this observation, we can formulate a spatial motion PDE as follows: 

Ut+Vx- [ua{u)] = (6) 

where a{u) is the velocity field of the spatial transport. As elaborated by the examples above, the spatial 
velocity field a{u) should be defined as 

a{u) = -^^{A^{Au - /)) • sign(u) 
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Thus ^ can be equivalently formulated as 



ut^V^-[\u\S/^{A^{Au^ f))]. (7) 

For definiteness we attach the initial and boundary conditions to ([7]) and obtain the following Cauchy problem 

ut^V,-[\u\V,{A^{Au~m, {x,t)enx{0,T] 
u{-,0) = uo,u\9n = 

where Uq is a continuous initial function. ([5]) is the core spatial motion PDE we will discuss throughout this 
paper. A locally integrable function u is called a weak solution of ([5]) if u satisfies the boundary condition 
and for any test function 77 the following identities hold: 

/ [7^fU-\u\V.,7^-V,iA^iAu- f))]dxdt^O, (9) 

JSlx(0,T] 

lim / rj(x,t)u{x,t)dx — I ri(x,0)uo(x)dx. (10) 
*->o+ Jn Jn 

3.2 Physical and Geometric Interpretations of PDE ([7]) 

The equation ([7]) has a nice physical interpretation, and is closely related to optimal transportation and 
nonlinear dissipative processes. To see this, we assume u > for now and regard u as the mass density of a 
certain fluid. The mass conservation can be expressed as the following continuity equation for u: 

Ut + \7^- [uw] = (11) 

rn 

where w is the velocity field that describes the spatial motion of the fluid. Darcy's law ^41] connects the 
velocity field w with the pressure field p by 

w = -V^p, (12) 
and thermodynamics tells us the pressure p is determined by the potential E{u): 

p^duE{u). (13) 

For an ideal gas in a homogeneous porous medium, the potential E[u) is given by the free energy J Cu'^ . 



Thus equations ([TT|). ([T2|) and ([TB]) lead to the classical porous medium equation 4l|. If we replace the 
potential E{u) by our residual H{u), the combination of (|lip . and ([T^ turns out to be our spatial 
motion equation ([7|). 

Equation JT]) can also be understood as a gradient flow of the potential H{u) under the Wasserstein 



metric. In |42l . l43j . Otto et al. showed that the porous medium equation is a gradient flow under the 
Wasserstein metric for some given potential. Our following interpretation of ([7]) will be in the same spirit 



as m 



42 



43|. 



Wasserstein distance is widely used in optimal transportation problems, whose fundamental goal is to 
find the most efficient plan to transport a density function to another density function ui in a mass 
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preserving fashion. Monge's optimal transportation problem is to find the solution to 



inf / uo{x)\x ~ M{x)\^dx (f4) 



44 



451 for the definition 



and the corresponding optimal value is called Wasserstein distance of order 2 (see 
of and more other details). The idea of optimal transportation has been applied to various problems in 



image and signal processin g, su ch as image classification [46|, registration 47| and segmentation [48| . 

As pointed out by Otto [42] , the Wasserstein distance of order 2 can be understood as a geodesic distance 
under a certain Riemannian structure on the set of density functions. The metric tensor g on the tangent 
space at u is formally defined by 



(si,S2)g= J U^x4'l ■^x(t)2 (15) 

where Si and S2 are any two infinitesimal variations of u, and i — 1, 2, are solutions of 

S^+V^■ {U\/^(|)^) ^0. (16) 

Given the metric tensor g, the set of all density functions forms a Riemannian manifold {A4,g). 

Now we can formally interpret the motion PDE ([7]) as a gradient flow of H{u) = — /HI on {A4,g). 

Under the Riemannian structure, the corresponding gradient of energy H{u) w.r.t. g, denoted by gia.dgH{u), 
is defined by the following identity 

{gradgH{u),v)^^d,^H{u) (17) 

for all vector field v on {A4,g), where d^j denotes the directional derivative along v. On the other hand, from 
([Tsl) and (fT6|) we can sec that (si, 52)^ = / 0iS2 where ipi solves si + Vx ■ {u'Vx4'i) ~ 0- Therefore, since u 
solves ut + Vx ■ {u'Vx{—duH{u))) = 0, we have 



{ut,v)g = J -duH{u) -v^-d^Hiu), W (18) 

which indicates that Ut — — gradgi?(it). This elegant interpretation leads to a substantial understanding of 
our motion PDE ([7]) : it gives the most natural way (in the sense of optimal transport) to spatially move the 
signal such that the residual H{u) is decreasing. 



3.3 Properties 

Now we come back to the general case where u is not assumed to be positive. We have the following fairly 
well known properties of equations that resemble ([S]) that explain, to some extent, why this PDE can be 
utilized to improve the spatial reconstruction. 

Proposition 1. Ifu{x,t), {x,t) G il x {0,T] is a weak solution of (|S]) then H{u{x,t)) = ^\\Au — /Hj is 
non-increasing over time. 
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Proof. We need to show that —H{u{x,t)) < 0. Indeed, 



dt 



H{u{x,t)) — / H' {u) ■ utdx 



V^iA'iAu - /)) • [\u\VM^{Au - f))]dx 
J \u\[V^iA'^iAu- f))]^dx<0 



□ 



Proposition 2. Ifu{x,t), {x,t) gQx (0,T] is continuous weak solution of ([8]) then \\u\\i is conserved over 
time. 

Proof. For simphcity, we will prove this in the ID case, while the proof for higher dimensions is similar. For 
a fixed time t, suppose / = [a, fe] is one of the connected component of {u ^ 0}, i.e. u has constant sign 
within /, then we have u(a, t) = u{b, t) = 0, so 



d 
di 



\u\dx 



utsigii{u)dx 



nb 

signH- / V, •[|u|V,(^^(A«-/))] 

J a 

sign{u)-[\u\V4A'^{Au-f))]\lZl = 



Therefore overall we have = 



□ 



The above two propositions show that, although the total mass of u is preserved over time, the spatial 
mass distribution is adjusted such that the equation Au = / is better satisfied. 



o 


□ 










o 


o 


□ 










c 






o 






El 



Figure 1: The first column are observed images / obtained from convolution of a Gaussian with one or 
two spikes. The the second column shows the initial functions uq and the remaining 4 columns show the 
corresponding solutions of spatial motion PDE (O at time T=250, 500, 750, 2000 respectively. The difference 
between the second and third rows reflect the fact that with the same /, different choices of initial value will 
lead to different results. 
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Another and yet more important property of the motion PDE is that the sohition u{x,t) becomes 
sparser asymptotically in t. We illustrate this phenomena numerically in Figure [l] where the peaks of the 
solution u move to the correct locations, and when they are close to the desired locations, u begins to become 
sparser which is desired for sparse reconstruction problems. This explains why we can get more accurate 
solutions when combining the motion PDE (|S|) with FPC Bregman iterations (see Section 14.3.21 for more 
details). 

It is also worth noticing that, comparing the second and the third rows of Figure [TJ for different initial 
functions uq, the solutions of PDE (|8]) have rather different asymptotic behavior. This means that the 
residual H{u) is not convex under the Wasserstein metric. In fact, the solution in the second row of Figure 
[T]only approximates a local minimizer of H{u) with respect to the Wasserstein metric but not a global one, 
i.e. H{u) ^ in this case, but gradgif(u) — 0. 



4 Numerical Schemes and Experiments 



4.1 Numerical Implementation of PDE ([H]) 

The finite difference scheme that we use for our motion PDE is as follows: 



,n+l 



dt 



Dlpl 



(19) 



where 









,j — Ui+lJ Uij, 






Dlu 




( K 


2 




Ui+ijUij > 0, 


"I 






otherwise 


r K 


2 




Uij+iUij > 0, 


"I 






otherwise 



and 



where the operator A is implemented by standard matrix multiplication. The special definition of ^ | 
and I above is to make sure that the £i norm of u'^ is exactly preserved in the discrete setting. 

To guarantee the stability of the numerical solution, the time step dt should be chosen appropriately. 
Although it is hard to derive a explicit stable condition for our nonlinear PDE, for each step we can treat 
the PDE as a transport equation with a fixed velocity field, and find a necessary stability condition, which 
is given by 



1 



>me.x{\Dlpl^\,\Dlp-J}. 



(20) 



4.2 Algorithm 



As we explained in section l373l the PDE ([8]) does not necessarily lead to a correct solution of the optimization 
problem ([1]). However, it helps the function to be more sparse while decreasing the residual. Therefore it can 
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be combined with other existing deconvolution algorithm and improve the performance. Since the motion 
PDE ([HI is time-dependent, it is natural to combine the numerical scheme p9|) with FPC Bregman iteration 
(HI) and introduce the following algorithm: 



u' 



.k+l,N 




yfcj+i ^ shrink(w''''+t,^) 



(PDE Update) 



(21) 



fk+l ^fk^f_ ^ufc+i,JV_ 

We shall refer to the above algorithm as PDE-FPC Bregman iterations. 

We further note that for each iteration the PDE update only slightly increases the complexity, because we 
do not need to recalculate the term A^ {Au — /) and the finite difference operator ([T9| can be implemented 
rather efficiently. 

We want to point out that this PDE update step can be similarly imposed to many other iterative 
thresholding based algorithms, e.g. the linearized Bregman iterative method [2^ or the BOS method [s^. 
As far as we tested, this complementary step always improves the performance and accuracy of these algo- 
rithms without increasing the computational complexity very much. In many cases the improvement is very 
significant. 

Remark: Note that the inner loop of (I^Tt contains TV > 1 iterations, as in the original FPC algorithm. 
The choice of N is not critical for the convergence of the scheme, however it affects the performance. When 
N is small, the algorithm converges faster (in terms of total number of iterations) while the decay of the 
residual is more oscillatory; when N is large, it converges slower while the decay of the residual appears more 
stable. In our experiments we choose = 10, which seems to give us a good balance between the speed and 
decay of residual. 

4.3 Numerical Results 

In this section, we compare the performance of FPC Bregman interations (HJ with that of PDE-FPC Bregman 
iterations (j2ip . We will first consider the case where there is not any noise in /, and compare the convergence 
speed. Then we consider a noisy case and compare the accuracy of the solutions obtained from the two 
algorithms. The operator A is taken to be a convolution operator with a Gaussian kernel. The clean and 
noisy / are shown in Figure [2l where the size of the images is 50 x 50. The two kernels of the clean 
blurry images are generated by Matlab function f specialC 'gaussian' , [41 41] ,a) with cr = 4 and 4.5 
respectively. The noisy images are generated by adding Gaussian white noise to the clean blurry image with 
kernel f speciaK 'gaussian' , [21 21] ,5.5). The intensities of the spikes are generated randomly in [0.5, 1] 
with locations indicated by the black circles (Figured]). 

4.3.1 Convergence Speed 

We consider the noise free case (left two images in Figure [2]). Here we take the stopping criteria < 
10~^, where u is the true solution. The following table summaries the comparison results with various 
choices of parameter fj,. The time step dt for the PDE is chosen to reflect best performance for given fj,. 
From Table [T] and Table [H we can see that PDE-FPC Bregman is generally faster than FPC Bregman in 
terms of computational time. Furthermore, the smaller ^ gets, the more PDE-FPC Bregman outperforms 
FPC Bregman. The reason is that for small /i, u*'' is more regular (i.e. less sparse) than for large fj, due to 
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Figure 2: Clean observed images generated by different Gaussian kernels (left and middle). Noisy observed 
image (right) with SNR=15.87. Black circles indicating the locations of the spikes. 

thresholding, and hence the PDE update in (PT|) will have a much nicer behavior. Also, a larger time step 
dt is allowed when u'^ is more regular. 

It is also worth noticing from Table [T] and [2] that the best performance of PDE-FPC Bregman is 6 and 
10 times faster respectively than that of FPC Bregman (marked by "♦"). The more blurry is the image, the 
bigger is the improvement of PDE-FPC Bregman over FPC Bregman. 

Table 1: Comparisons of FPC Bregman Q and PDE-FPC Bregman The symbol "4" labels the best 

results for FPC Bregman and PDE-FPC Bregman. 







FPC Bregman 


PDE-FPC Bregman 




dt 


^Alterations 


Time (sc) 


^Iterations 


Time (sc) 


(2, 0.5) 


0.2 


2824 


13.21 


2816 


14.97 


(2, 0.2) 


0.5 


2196 


10.69 4 


1705 


9.23 


(2, 0.1) 


2.0 


2585 


12.43 


1324 


7.06 


(2, 0.05) 


4.5 


3365 


16.04 


1291 


6.91 


(2, 0.01) 


300.0 


4997 


24.20 


587 


3.20 


(2, 0.002) 


800.0 


9622 


47.76 


280 


1.64 


(2, 0.001) 


1200.0 


13386 


67.73 


243 


1.48 4 



4.3.2 Accuracy 

We now consider the noisy case (right image in Figure [5]). Experimental results are summarized in the 
following Table [31 Since there is noise in /, the error, which is taken to be ^^jh]|7^' generally decreases 
first and then increases, and will eventually converge to something noisy. Therefore, in order to reflect 
the best performance of both algorithms, we recorded the number of iterations and computation time that 
correspond to the smallest error possible for each given set of parameters. From Table [3] we can see that, 
for each given /i, although PDE-fFPC Bregman may not always be faster than FPC Bregman, it always 
produces a solution with a smaller error. If we consider the best set of parameters for the two algorithms, 
PDF-FPC Bregman outperforms FPC Bregman in terms of both computation time and accuracy. 
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Table 2: Comparisons of FPC Bregman (g]) and PDE-FPC Bregman (HJ). The symbol labels the best 
results for FPC Bregman and PDE-FPC Bregman. 







FPC Bregman 


PDE-FPC Bregman 




dt 


^Iterations 


Time (sc) 


^Iterations 


Time (sc) 


(2, 0.5) 


0.25 


5327 


25.37 


4808 


25.33 


(2, 0.2) 


0.9 


4115 


19.81 4 


2904 


15.26 


(2, 0.1) 


2.0 


4651 


22.13 


2562 


13.70 


(2, 0.05) 


5.0 


5836 


28.31 


1876 


10.31 


(2, 0.01) 


300.0 


8460 


41.42 


1030 


5.61 


(2, 0.002) 


1500.0 


14996 


74.01 


446 


2.53 


(2, 0.001) 


2000.0 


22466 


112.29 


328 


1.91 4 



Table 3: Comparisons of FPC Bregman (g]) and PDE-FPC Bregman The symbol labels the best 

computational time and error respectively for FPC Bregman and PDE-FPC Bregman. 







FPC Bregman 


PDE-FPC Bregman 




dt 


^Iterations 


Time (sc) 


Error 


^Alterations 


Time (sc) 


Error 


(2, 0.8) 


0.01 


13299 


17.52 


1.12e-001 


13417 


24.61 


1.04e-001 


(2, 0.5) 


0.05 


6086 


7.74 


2.68e-002 4|k 


6033 


10.81 


1.63e-002 


(2, 0.2) 


0.5 


3663 


4.61 


5.70e-002 


3321 


5.94 


2.40e-002 


(2, 0.05) 


1.5 


2127 


2.94 <|k 


4.87e-002 


1424 


2.59 


1.42C-002 


(2, 0.01) 


3.0 


2404 


3.49 


5.17e-002 


1651 


3.09 


1.47e-002 


(2, 0.005) 


4.0 


2511 


3.81 


9.00e-002 


1406 


2.73 


1.93e-002 


(2, 0.002) 


45.0 


2985 


4.92 


1.87e-001 


571 


1.20 4 


1.26C-002 4 


(2, 0.001) 


40.0 


2730 


4.70 


4.05e-001 


623 


1.38 


5.22e-002 



5 Conclusions 

In this paper, we introduced a new nonlinear evolution PDF for sparse deconvolution problems. We showed 
that this spatial motion PDF preserves the £i norm while lowering the residual \\Au — /||2- We also showed 
numerically that the solution of the motion PDF tends to become sparse and converges to delta functions 
asymptotically. Therefore, utilizing these properties of the PDF, our proposed PDE-FPC Bregman itera- 
tions (|2T|) outperformed the original FPC Bregman iterations [l| in terms of both convergence speed and 
reconstruction quality. This was strongly supported by our numerical experiments. We finally note that our 
PDE can also be used to enhance performance for many other iterative thresholding based algorithms, e.g. 

3- 



25 



2a 
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